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Abstract 

This  thesis  presents  an  algorithm  suitable  for  numerical  analysis  of  cryogenic  refrigeration 
systems. -The  need  for  such  an  algorithm  arises  from  the  absence  of  effective  analysis  tools  in 
the  cryogenic  field. 

Typically^- preliminary  design  of  a  cryogenic  system  commences  with  a  number  of  decoupling 
assumptions  with  regard  to  the  process  variables  of  heat  and  work  transfer  (e.g.  work  input  rate, 
heat  loading  rates)  and  state  variables  (pinch  points,  momentum  losses).  These  assumptions  are 
made  to  facilitate  direct  and  simplified  solution  calculations.  However, "making  preliminary 
performance  estimations  minimizes  the  effect  of  component  interactions  which  is  inconsistent 
with  the  intent  of  "analysis" .  A  more  useful  design  and  analysis  tool  is  one  in  which  no  restriction  s 
are  applied  to  the  system  -  interactions  become  absolutely  coupled  and  governed  by  the  equi¬ 
librium  state  variables.  Such  a  model  would  require  consideration  of  hardware  specifications 
and  performance  data  and  information  with  respect  to  the  thermal  environment.  Model  output 
would  consist  of  the  independent  thermodynamic  state  variables  from  which  process  variables 
and  performance  parameters  may  be  computed.  This  model  will  have  a  framework  compatible 
for  numerical  solution  on  a  digital  computer  so  that  it  may  be  interfaced  with  graphic  symbology 
for  user  interaction. 

This  algorithm  approaches  cryogenic  problems  in  a  highly-coupled  state -dependent  manner. 
The  framework  for  this  algorithm  revolves  around  the  revolutionary  thermodynamic  solution 
technique  for  Computer  Aided  Thermodynamics  (CAT),  developed  by  Dr.  Gilberto  Russo  at 
the  Massachusetts  Institute  of  Technology. 'Fundamental  differences  exist  between  the  Control 
Volume  (CV)  algorithm  and  CAT,  which  will  be  discussed  where  appropriate. 

This  thesis  presents  the  algorithm  in  detail  with  respect  to  the  thermodynamic  principles, 
state-variable  management,  device  models,  system  networking,  solution  structure,  numerical 
reduction,  model  flexibility  and  generality  and  potential  pitfalls. 
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1  Introduction 


1.1  Background 

Much  design  work  has  been  done  in  the  cryogenic  engineering  field  and,  accord¬ 
ingly,  there  are  numerous  design  methods  that  have  proven  useful.  Design,  in  an  engi¬ 
neering  context,  means  to  model  physical  behavior  and  assign  value  to  problem  variables 
in  such  a  manner  as  to  achieve  a  desired  engineering  effect.  A  "good"  design 
distinguishes  itself  from  a  "bad"  design  by  whether  or  not  the  variables  and  parameters 
were  assigned  value  in  an  optimal  sense.  Optimal  is  relative,  but  usually  implies  maxi¬ 
mum  desired  effect  at  minimal  resource  cost.  Efficiencies  and  dimensional  analysis  vari¬ 
ables  are  typical  ways  in  which  engineers  evaluate  the  optimization  of  a  particular  design. 
Thus,  good  design  requires  qualify  optimization.  Quality  optimization  requires  a  thorough 
understanding  of  the  problem  variables,  or,  more  succinctly,  good  analysis. 

Analysis  means  understanding  the  interactive  behavior  of  the  relevant  physical  vari¬ 
ables.  This  is  the  reason  good  design  work  is  supported  by  a  sound  foundation  of  analyti¬ 
cal  study.  Good  analysis  tools  are  not  in  abundance  in  the  cryogenic  field,  which  forms 
the  basis  for  this  study:  to  develop  an  analysis  algorithm  suitable  for  application  to 
cryogenic  engineering  systems.  A  typical  analysis  method  for  a  closed  loop  thermody¬ 
namic  system  starts  with  a  definition  of  the  physical  system  in  terms  of  models  for 
devices  (turbines,  boilers,  etc.)  and  for  their  interaction  with  the  environment  (Q,n,  Wout, 
TRHS,etc.) .  To  determine  the  system  performance,  a  number  of  assumptions  must  be  made 
to  reduce  the  vast  number  of  unknowns.  Such  assumptions  include  the  thermodynamic 
state  at  the  at  the  inlet/outlet  of  a  device  (e.g.  saturated,  superheated,  quality),  a  ratio  of 
relevant  pressures  (AP/P,  PH/P, ),  or  a  temperature  change  ratio  (T„/T[ ,  pinch  points).  The 
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effect  of  making  these  assumptions  is  to  decouple  or  minimize  the  interaction  of  state 
variables  that  govern  system  behavior,  the  behavioral  patterns  of  the  system  become 
locked  out. 

A  robust  analysis  is  one  which  minimizes  global  assumptions  allowing  the  process 
variables  to  result  directly  from  the  highly  coupled  interaction  of  equilibrium  state  vari¬ 
ables.  The  analysis  method  could  be  likened  to  a  black  box  which  takes  inputs  of  compo¬ 
nent  specifications  and  yields  as  output  the  value  of  the  system's  equilibrium  state 
variables.  It  is  then  a  simple  extension  to  transform  the  state  variables  into  more 
significant  performance  parameters.  Thus  system  modelling  is  achieved  by  allowing  the 
system  to  operate  in  its  unrestricted  condition  compatible  with  the  hardware  imposed. 
Optimization  (and  hence  design)  becomes  a  matter  of  selecting  hardware  that  achieves 
the  best  overall  system  performance,  as  defined  by  the  designer.  The  method  which  is  the 
heart  of  this  paper  pursues  this  strategy. 

Analysis  is  difficult.  It  requires  simultaneous  segregation  and  synthesis  of  the  prob¬ 
lem  variables  to  establish  process  characteristics  or  system  behavior.  Many  engineers 
incorrectly  use  the  word  analysis  when  referring  to  design  calculations.  In  fact,  analysis 
in  its  crudest  (yet  often  practiced)  form  is  executed  by  conducting  a  series  of  design  cal¬ 
culations  with  differing  input  values  and  observing  the  output  or  performance  values. 

This  is  analogous  to  the  experimenter  plotting  the  results  of  a  run  without  first  calculating 
the  analytical  prediction.  Such  an  analysis  method  does  the  segregation  while  loosing 
sight  of  the  synthesis.  This  type  of  investigation  is  more  appropriately  termed  design 
analysis. 
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1.2  Basis  for  Method 

In  this  work,  the  methodology  for  analysis  of  a  thermodynamic  system  consists  of 
idealization  of  the  physical  problem,  generation  of  the  mathematics  modelling  the  prob¬ 
lem,  solution  of  the  mathematical  representation  and  interpretation  of  the  results.  The 
aggregate  of  these  tasks  is  rather  complex.  The  method  presented  here  reduces  the 
complexity  by  using  simple  concepts  to  construct  a  complex  representation  of  thermody¬ 
namic  behavior  that  will  be  managed  by  a  computer.  This  is  essentially  what  is  done  in  a 
finite  element  analysis  method:  building  a  complex  representation  to  be  managed  by  a 
computer  from  almost  trivially  simple  engine. ring  concepts. 

The  great  utility  of  a  digital  computer  is  that  it  can  be  programmed  to  perform  the 
equation(s)  generation  and  mathematical  solution  (through  an  intelligent  assembly  cf  the 
"simple  pans"),  leaving  the  idealization  and  interpretation  to  the  use;  (although  some  of 
today’s  modem  routines  will  assist  in  the  interpretation  process  by  a  convenient  assembly 
of  the  output  data).  The  method  to  be  developed  will  allow  a  graphic,  symbolic  input  of 
the  problem,  thus  funher  assisting  the  user  in  the  analysis  of  alternatives.  The  ability  of  a 
computer  to  manage  large  quantities  of  data  (specifically,  unknown  variables)  eliminates 
the  temptation  to  make  decoupling  assumptions  about  problem  variables  that  reduces 
their  number  to  an  (algebraically)  manageable  level.  It  is  this  special  feature  that  makes 
this  method  unique:  elimination  of  decoupling  assumptions  with  respect  to  the  problem 
variables  causes  the  variables  to  become  mathematically  unconstrained  and  therefore 
physically  interactive  through  the  interconnective  matching  conditions. 

In  applying  a  computer  to  perform  thermodynamic  analysis,  the  method  employed 
should  be  general  enough  to  allow  analysis  of  a  wide  variety  of  complex  multiple-element 
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thermodynamic  systems.  In  such  an  employment,  it  should  use  simple  models  (vis.  equa¬ 
tion  sets)  to  represent  the  composite  system  and  the  same  solution  procedure  to  solve  for 
the  requisite  data.  This  was  the  inspiration  behind  Computer  Aided  Thermodynamics 
(CAT)  as  developed  by  Dr.  Gilberto  Russo  at  the  Massachusetts  Institute  of  Technology. 
CAT  is  structurally  similar  to  a  finite  element  analysis  except  that  it  is  applied  to  a  collec¬ 
tion  of  discrete  elements  (that  may  undergo  different  types  of  interactions)  rather  than  a 
homogeneous,  continuous  system  (undergoing  a  single  interaction).  The  great  success  of 
CAT  has  served  as  a  backdrop  for  the  Control  Volume  (CV)  algorithm:  adapting  thermo¬ 
dynamic  analysis  to  take  advantage  of  computer  power  by  formulating  a  new  method  of 
analysis. 


2  Foundations  of  Control  Volume  Method 

2.1  Overview 

The  control  volume  analysis  method  is  developed  for  implementation  on  a  digital 
computer.  The  basis  for  the  method  is  modelling  of  complex  systems  using  simple  build¬ 
ing  elements.  These  elements  come  in  two  broad  classifications:  storage  elements  and 
interconnective  elements.  Each  element  has  associated  with  it  a  set  of  constitutive 
relations  that  describe  the  element’s  behavior  and  are  defined  in  a  problem  independent 
manner. 

The  complete  physical  system  is  represented  by  all  the  constitutive  relations  of  all 
the  elements  making  up  the  system  specified  for  the  particular  arrangement  of  elements. 
Once  the  elemental  constitutive  relations  become  specified  for  the  system  (through  a 
change  of  base  or  variable  transformation),  they  become  the  basis  for  the  residual  rela¬ 
tions  of  the  system.  Residual  relations  are  used  to  define  the  steady  state  behavior  of  the 
entire  system.  The  specificity  comes  from  invocation  of  the  system  topology.  The  system 
topology  is  used  to  perform  the  change  of  base  that  generates  the  problem  specificity 
from  the  "n"  sets  of  constitutive  relations  of  the  "n"  elements  of  the  mesh.  This  accom¬ 
plishes  the  ’modelling  and  equation  generation’  steps  in  this  computer-aided  thermody¬ 
namic  analysis.  Solution  is  found  by  using  the  residual  functions  in  a  global  equilibrium 
equation: 

[K'(x,n{A.x}  ={R(x,)}  (1) 

(K'(x)l  is  a  matrix  of  first  derivatives  of  the  residual  relations  {R(x,))with  respect  to  the 
independent  variables  { x, ) of  the  system  (i.e.  the  Jacobian  matrix).  Successive  modifica- 
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tions  are  made  by  imcrementing  the  current  value  of  { x }  by  {Ax}  which  will  drive 
[K'fXjlHAx)  to  zero  (or  less  than  some  appropriate  e).  [Kl(x,)J{Ax}=0  satisfies  the 
constitutive  relations  of  the  system.  The  independent  variable  vector  { x, }  for  which 
{ R(x) }  goes  numerically  close  to  zero  is  the  equilibrium  state  of  the  system  and  no  fur¬ 
ther  increment  in  x  is  required. 

[K‘(x,)J  is  defined  as  the  System  Tangent  Stiffness  Matrix1.  When  (Kl(x,)]  is 
assigned  a  numerical  value  using  {x,},  [K'(x,)]  is  the  linear  approximation  to  the  system's 
residual  functions  {R(x,)}.  This  linearization  is  necessary  due  to  the  non-linear  nature  of 
some  of  the  system’s  equations.  The  System  Tangent  Stiffness  Matrix  represents  an  (n-1) 
dimensional  tangent  surface  to  a  (n)  dimensional  function,  hence  the  tag  "tangent".  The 
modifier  stiffness  comes  from  the  congruence  between  equation  (1)  and  the  equilibrium 
equation  for  a  structural  finite  element  analysis  system  where  the  constitutive  matrix  is 
often  referred  to  as  the  "stiffness"  matrix. 

This  paper  describes  the  physical  idealization  of  the  thermodynamic  problem,  gen¬ 
eration  of  the  mathematics  describing  the  problem  and  solution  procedure.  Interpretation 
of  results  is  not  considered  here,  since  that  is  a  matter  of  recasting  equilibrium  state 
variables  into  more  pragmatic  (or  new)  performance  representations  (COP.  r|m,  etc.)  that 
may  be  defined  in  the  post-processing  routines. 


1  "A  New  Methodology  of  Computer-Aided  Thermodynamics",  Gilberto  Russo  [doctoral 
thesis],  Department  of  Nuclear  Engineering,  Massachusetts  Institute  of  Technology,  January. 
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2.2  Compatibility 

The  algebraic  uniqueness  and  simplicity  of  each  storage  element  lead  one  to  believe 
that  system  simulation  is  now  just  a  matter  of  connecting  together  elements  of  interest 
and  performing  the  system  integration  and  reduction.  However,  the  issue  of  compatibility 
must  be  addressed  with  respect  to  the  interaction  of  the  engineering  devices  (  which  are 
independently  defined).  In  other  words,  does  the  ’whole  equal  ihe  sum  of  the  parts'. 

This  problem  commonly  arises  during  the  design  of  thermal  power  systems  with 
regard  to  "matching"  of  components.  For  example,  does  the  mass-flow/pressure  ratio 
characterisdc  of  a  compressor  match  that  of  the  turbine  it  is  to  operate  with?  If  not,  it  may 
be  impossible  to  achieve  a  steady  state  operation  for  the  system  in  the  absence  of  spe¬ 
cially  designed  controls.  Such  a  system  (without  controls)  may  oscillate  (its  performance 
and  state  variables)  about  a  stable  equilibrium  (i.e.  hunting),  one  component  may  drive 
the  other  to  operate  at  an  inefficient  operating  point.  The  stability  of  the  operating  state 
will  require  further  consideration  of  the  dynamics  of  the  system. 

Accordingly,  this  methodology  investigates  the  steady  state  operation  of  systems  for 
which  performance  parameters  (constants  or  functions)  are  assigned.  Furthermore,  the 
relations  and  models  used  to  describe  component  behavior  must  be  consistent  with  physi¬ 
cal  limitations  of  a  system.  For  example,  it  is  incompatible  to  specify  pressure  rises  for 
both  a  compressor  and  an  expander  which  operate  in  the  same  fluid  circuit  since  this 
would  be  to  over-constrain  the  problem.  Only  one  pressure  ratio  can  be  specified  while 
the  other  pressure  ratio  results  as  a  consequence  of  the  interaction  of  the  pressure  and 
mass  flow  variables  in  the  system.  This  compatibility  is  necessary  to  ensure  a  steady  state 
can  be  attained. 
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To  accommodate  the  compatibility  issue  during  the  initial  stage  of  the  development 
of  the  method,  pressure  ratios  will  be  defined  only  for  compression  devices.  All  expan¬ 
sion  devices  (turbines,  throttles,  etc.)  will  be  void  of  any  pressure  characteristic;  the  equi¬ 
librium  pressures  will  be  established  by  the  interaction  of  the  system  elements. 

2.3  Physical  Idealization 

The  physical  system  is  modelled  using  its  constitutive  elements  .  Figure  1  is  a  repre¬ 
sentation  of  a  closed  loop  thermodynamic  circuit  operating  as  a  refrigeration  cycle.  While 
heat  and  work  interactions  are  an  intimate  part  of  the  complete  operating  system,  these 
external  interactions  are  determined  after  the  equilibrium  behavior  of  the  system  is  estab¬ 
lished.  The  work  or  heat  interaction  for  any  element  is  readily  determined  by  simply 
knowing  the  inlet  and  outlet  state.  This  is  infact  the  approach  used  here:  work  or  heat 
transfers  are  determined  after  the  value  of  the  equilibrium  state  variables  are  known  by 
writing  the  1st  Law  of  thermodynamics  in  the  post  processing  routines. 

The  system  is  the  collection  of  control  volumes  that  contains  all  components  of 
interest  and  the  environment  with  which  the  components  interact.  A  valid  thermal  analy¬ 
sis  of  the  component  interactions  can  only  be  done  by  analyzing  a  uniquely  independent 
set  of  control  volumes.  This  system  has  four  distinct  components2.  Accordingly,  a  set  of 
four  independent  control  volumes  must  be  developed.  A  selection  of  any  three  of  the  four 
control  volumes  shown  in  Figure  1  plus  a  control  volume  encompassing  the  entire  system 
(including  the  environment,  which  is  part  of  the  system)  is  a  valid  selection.  However, 


2  For  simplicity  in  this  development,  interconnective  piping  and  flow  passages  are  assumed 
to  be  lossless.  Lossy  flow  in  interconnecting  piping  may  be  incorporated  into  the  framework 
upon  extension  of  the  applicable  constitutive  relations. 
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the  control  volume  methodology  is  based  on  using  simple  construction  elements  to  build 
the  total  system  model,  therefore  the  four  component  control  volumes  (including  their 
external  interactions)  idealize  the  system. 

In  choosing  the  relations  that  will  be  sufficient  to  describe  the  system  being  ideal¬ 
ized,  careful  thought  must  be  employed.  If  an  engineering  system  is  not  completely  speci¬ 
fied  in  terms  of  the  relevant  behavioral  relations,  then  the  number  of  additional  arbitrary 
parameters  that  must  be  specified  will  equal  the  number  of  degrees  of  uncertainty  left  in 
the  mathematical  representation  of  the  physical  system.  On  the  other  hand,  if  the  system 
is  overspecified,  ambiguous  and  contradictory  results  will  precipitate  (i.e.  over¬ 
constrained,  physical  incompatibility).  There  must  be  a  compatibility  between  the  physi¬ 
cal  and  mathematical  problem:  the  number  of  independent  unknowns  must  equal  the 
number  of  independent  relations.  One  of  the  obstacles  of  this  effort  was  in  finding  this 
delicate  balance  between  the  mathematical  and  physical  requirements. 
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2.3.1  Local  and  Global  Variables 

Local  variables  are  independently  defined  within  each  storage  element  and  are 
used  in  the  constitutive  relations  of  each  storage  element.  When  elements  are 
assembled  together  to  form  a  system,  the  local  variables  are  now  transformed  into  a 
variable  system  defined  for  the  complete  physical  system.  This  system  of  variables  is 
called  the  global  variables  system  and  it  is  defined  only  for  a  specific  system  based  on 
the  specific  interconnections  between  storage  elements. 

2.3.2  Storage  Elements 

A  storage  element  is  in  a  pragmatic  sense  a  self-contained  standard  engineering 
component  such  as  a  turbine,  heat  exchanger  or  valve  that  is  defined  on  an  elemental 
basis.  Storage  elements  transfer  energy,  entropy,  and  mass  to  achieve  some  desired 
engineering  effect.  Accordingly,  the  1st  Law  of  Thermodynamics,  continuity  require¬ 
ments,  and  an  efficiency  (irreversibility  parameter)  relation  are  chosen  to  describe  the 
constitutive  behavior  of  the  element. 

For  devices  which  exchange  work  with  their  environment,  the  classical  defini¬ 
tions  of  component  isentropic  efficiency  serves  as  the  irreversibilty  parameter.  Such  a 
definition  reflects  2nd  Law  irreversibilities  while  incorporating  component  specific 
"test-stand"  data. 

The  irreversibility  relation  for  devices  which  exchange  heat  (either  internally  or 
with  the  environment)  is  from  the  heat  transfer  rate  equation.  The  irreversibility  infor¬ 
mation  is  contained  in  the  overall  heat  transfer  coefficient  (U)  and  the  available  heat 
transfer  area  (A)  which  are  used  in  the  rate  equation  to  relate  fluid  stream 
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temperatures.  A  second  irreversibility  parameter  is  from  the  fluid  friction  losses.  Thus 
relations  defining  the  pressure  drop  through  a  heat  exchanger  together  with  the  heat 
transfer  rate  equation  define  the  irreversibility  (entropy  generation).  This  model  may 
seem  to  make  the  methodology  more  complex,  but  it  achieves  the  goal  of  incorporat¬ 
ing  hardware  performance  data  into  the  elemental  models. 

The  equations  used  to  specify  the  irreversibilties  are  called  the  characteristic 
relations  since  they  reflect  component’s  characteristic  behavior. 

The  working  fluid  constitutes  an  additional  constitutive  relation  for  each  fluid 
state.  The  type  of  working  fluid  employed  should  be  arbitrary  in  a  useful  analysis 
method.  The  physical  properties  of  the  fluid  should  be  expressible  as  a  function  of 
independent  variables.  This  has  been  accomplished  for  the  actual  fluid  states,  but  for 
the  isentropic  outlet  states  (work  transfer  devices),  working  fluid  properties  are 
expressed  as  a  function  of  the  independent  pressure  variable  and  the  value  of  the  isen- 
trope.  That  is,  houl=h(Poul,  sj. 

If  a  system  uses  a  single  working  fluid,  then  all  states  are  calculated  using  the 
same  property  relation.  On  the  other  hand,  if  circuit  contains  multiple,  independent 
fluid  loops  each  employing  different  fluids(e.g.  a  nitrogen  pre-cooler  for  a  liquid 
helium  plant),  then  a  fluid  designation  must  be  associated  with  each  loop.  In  other 
words,  each  fluid  loop  would  have  a  working  fluid  constitutive  relation  applicable  to 
that  fluid.  This  is  more  efficient  than  defining  the  fluid  constitutive  behavior  for  every 
element. 
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In  summary,  each  storage  element  has  constitutive  relations  grouped  categori¬ 
cally  by  1st  Law,  characteristic  relations,  mass  continuity  and  working  fluid  constitu¬ 
tive  behavior. 

2.3.3  Interconnective  Elements 

The  interconnective  elements  model  (through  their  constitutive  relations)  the 
linking  of  various  storage  elements  which  constitute  the  system.  As  individual  compo¬ 
nents  are  assembled  with  their  associated  elemental  constitut  /e  relations  and  local 
variables,  a  transformation  is  necessary  to  shift  from  the  local  variables  system  to  the 
global  variables  system.  Such  a  transformation  is  accomplished  by  the  interconnective 
element’s  constitutive  relations.  By  transforming  from  a  local  to  a  global  variables 
system,  topological  information  is  implicitly  contained  in  the  new  variables  system. 
The  interconnectivity  requirements  (i.  e.  the  constitutive  relations  of  the  interconnec¬ 
tions)  enforce  on  the  system  the  equality  of  independent  local  variables  at  a  nodal  con¬ 
nection  between  storage  elements.  Since  these  local  variables  are  equal,  they  assume 
the  global  variable  label  at  the  node.  This  ensures  component  interaction.  For 
example,  in  Figure  1,  the  outlet  pressure  from  the  compressor  must  equal  the  inlet 
pressure  to  the  high  temperature  heat  exchanger.  The  two  local  variables  are  equal  and 
are  assigned  the  global  pressure  variable.  Similar  operations  are  performed  on  all  the 
independent  variables  (Chapter  3).  When  the  transformation  is  complete,  the  number 
of  global  variables  equals  the  number  of  local  variables  less  the  number  of  intercon¬ 
nection  constitutive  relations. 

2.4  Equation  Generation 

The  constitutive  relations  of  the  system  are  expressed  as  the  { R(x,) }  vector  which  is 
derived  from  the  constitutive  relations  of  the  elements  and  the  topology  of  the  system. 
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Equation  generation  refers  in  this  context  to  the  generation  of  the  [K‘(x,)]  matrix  and  the 
{ R(Xj) }  vector.  These  quantities  are  generated  through  recursive  relations  which  draw  on 
topological  information  to  generate  the  appropriate  elements  of  both  [Kl(x)]  and  {R(x)}. 
The  methods  for  generating  these  relations  are  discussed  in  further  detail  in  Chapter  6. 

2.5  Solution 

The  solution  technique,  briefed  in  Section  2.1  is  a  Newton’s  method  tailored  for  thi 
application  and  will  be  reviewed  in  further  detail  in  Chapter  7. 
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3  Variables  and  Parameters  Management 
3.1  Selection  of  Independent  Variables 

Section  2.2  made  indirect  reference  to  the  properties  of  temperature  and  pressure  as 
the  independent  thermodynamic  variables  of  this  method.  In  fact,  a  significant  portion  of 
the  CV  algorithm  development  was  performed  with  temperature  and  pressure  as  the  inde¬ 
pendent  variables.  However,  temperature  and  pressure  are  intensive  thermodynamic  prop¬ 
erties.  Two  intensive  properties  can  fix  the  thermodynamic  state  in  a  homogeneous, 
single  phase  system  only.  In  generalizing  this  algorithm  to  an  operating  refrigeration 
cycle,  determination  of  properties  and  states  in  a  binary  or  two  phase  system  is  necessary, 
since  much  of  the  refrigeration  effect  results  from  the  latent  heat.  Therefore,  to  lend  gen¬ 
erality  to  the  algorithm,  pressure  and  specific  volume  were  selected  as  two  independent 
variables. 

For  shaft  work  machinery(turbines,  compressors,  etc),  three  equilibrium  states  are 
of  interest.  The  equilibrium  inlet  state,  the  isentropic  outlet  state  and  the  actual  outlet 
state.  The  isentropic  outlet  state  can  be  determined  as  a  function  of  the  inlet  entropy  value 
and  outlet  pressure  value:  h2s=h(P2,s1),  where  s1=s1(P,,v1). 

In  addition  to  pressure  and  specific  volume,  mass  flow  rate  has  also  been  selected  as 
an  independent  variable.  Mass  flux  scales  the  system’s  energy  and  entropy  interactions. 
Mass  flow  rate  is  an  essential  variable  when  devices  such  as  splitting  valves  and  mix¬ 
ing/recombination  valves  are  employed. 
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These  three  variables  (pressure,  specific  volume,  mass  flow  rate)  taken  together 
with  the  interconnectivity  relations  (between  devices)  ensure  full  thermodynamic  match¬ 
ing  at  the  "nodes"  between  storage  elements.  Thus,  a  complete  engineering  link  is  made 
between  two  devices  exchanging  fluid  energy. 

3.2  Dependent  Variables 

The  dependent  state  variables  are  all  other  thermodynamic  properties: 

h(P,v) 

hs(P,s(P,v)) 

s(P,v) 

T(P,v) 

Work  and  heat  transfer  are  also  system  unknowns,  although  not  variables  perse.  They  are 
determined  by  the  equilibrium  values  of  the  global  variables.  The  1st  Law  residual  func¬ 
tions  and  [Kl]  entries  for  elements  which  exchange  heat  or  work  with  the  environment  are 
not  incorporated  as  part  of  equation  (1)  framework.  Such  relations  are  assigned  to  a 
post-processing  routine.  In  doing  this,  no  specificity  of  the  [K‘(x)]{Ax}  =  {R(x)}  equilib¬ 
rium  equation  is  lost  since  elimination  of  a  1st  Law  relation  also  eliminates  either  a  work 
transfer  or  heat  transfer  unknown. 

Determination  of  the  equilibrium  state  variables  (x)  completely  defines  the  state  of 
the  system,  from  which  heat  and  work  interactions  can  be  determined. 
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3.3  System  Parameters 

System  parameters  are  user  defined  as  pan  of  the  algorithm  initialization  where  the 
specifications  of  component  data  are  imported  for  use  in  the  applicable  constitutive  rela¬ 
tions.  Such  parameters  include  component  efficiencies  Cn.^Hc)-  friction  factors  (f),  heat 
transfer  coefficients  (U),  heat  transfer  areas  (A),  and  expansion/compression  ratios  (rE,  rj. 
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4  Device  Models  for  Control  Volume  Analysis 

4.1  General 

A  set  of  constitutive  relations  for  'he  individual  control  volumes  has  been  devel¬ 
oped  which  fit  into  a  standardized  framework  based  on  the  engineering  principles  dis¬ 
cussed  in  Chapter  2.  The  constitutive  relations  are  grouped  categorically  into  1st  Law. 
continuity,  interconnectivity  and  characteristic  relationships  (characteristic  relations 
describe  the  non-ideal  behavior  of  the  selected  control  volume  in  view  of  the  Second  Law 
irreversibilities  that  degrade  the  system  performance  from  the  ideal  case).  The  categorical 
standardizafon  of  the  relations  for  different  types  of  engineering  devices  is  necessary  to 
facilitate  the  assembly  of  the  elements  of  [K!( { x } )]  (Ax)={R({x))}  via  numerical  meth¬ 
ods. 

4.2  Engineering  Devices 

The  following  devices  have  been  selected  for  implementation  in  the  control  volume 
analysis: 

Regenerative  Heat  Exchangers 

Ambient  Heat  Exchangers 

Compressors 

Expanders 

Expansion  Valves 

While  this  list  is  not  all-inclusive,  it  is  representative  of  typical  devices  found  in 
cryogenic  refrigeration  systems  and  thermal  power  systems  in  general.  More  complicated 
devices  such  as  Stirling  engines  are  not  covered  here  but  may  be  included  upon  extension 
of  the  constitutive  framework  discussed  in  Chapter  2. 
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4.3  Component  Constitutive  Relations 

Symbology  for  the  expander  is  shown  in  Figure  2. 

lst  Law:  rii7h7(P2,v2)  =  WF 

Characteristic  relations: 

(hi-h2)-r\t(hu-h2s)  =  0 
hx-hu  =  0 

5^  —  5,  —  0 

The  first  equation  is  the  definition  for  isentropic  efficiency.  The  second  equation  is  a  con¬ 
sequence  of  the  recursion  used  to  generate  the  efficiency  equation:  h2s  can  not  be  numer¬ 
ically  generated  without  generating  his.  This  is  because  the  recursive  relations  use 
topological  information  to  generate  the  appropriate  functions.  Topological  information 
can  discriminate  between  inlet  and  outlet  states  but  not  between  isentropic  and  actual 
states.  The  recursion  generates  an  isentropic  inlet  state  which  is  identical  to  the  actual 
inlet  state;  the  second  equation  merely  sets  the  two  states  equal.  The  third  equation 
equates  the  inlet  entropy  to  the  outlet  entropy  for  an  isentropic  process(which  must  be 
known  to  use  the  efficiency  definition).  is  used  in  a  fluid  constitutive  relation  to  deter¬ 
mine  h2s=h(P2,  Si,). 

continuity:  =  0 

working  fluid:  h ^  -  h(P2,s,J  =  0 
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The  working  fluid  constitutive  relation  defines  the  outlet  enthalpy  for  the  isentropic  pro¬ 
cess. 


The  nine  unknowns  in  these  equations  (two  mass  fluxes,  two  pressures,  two  specific 
volumes,  an  isentropic  enthalpy  and  entropy  and  work  rate)  are  represented  by  only  six 
equations.  Since  each  control  volume  is  a  basic  building  element  of  the  CV  algorithm, 
each  individual  control  volume  should  be  algebraically  independent.  Therefore  three 
additional  equations  are  necessary  -  these  equations  come  from  the  interconnectivity 
requirements  for  pressure,  specific  volume  and  mass  rate.  Each  (2-port)  control  volume 
will  have  associated  with  it  three  interconnectivity  relations  (from  an  adjacent  intercon¬ 
nective  element)  which  add  to  the  other  constitutive  relations  to  ensure  algebraic  unique¬ 
ness. 

The  applicable  constitutive  relations  can  also  be  developed  for  the  compressor 
engine,  which  is  schematically  shown  in  Figure  3. 

1st  Law: 

rhlhl(Pl,v])-m2h2(P2,v2)  =  Wc 

Characteristic  relations: 

thermal  characteristic: 

-h2)~(hu-h2t)  =  0 
h\~hu  =  Q 

**-*i=0 
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pressure  characteristic: 


mass  continuity: 


working  fluid: 


P\~rcp  2  =  0 


mx-m2  =  0 


h7s~h(P2,s2s)  =  0 


These  equations  are  similar  to  those  for  the  expander,  accept  that  a  pressure  rise  relation 
(which  sets  the  system  pressure  levels)  has  been  incorporated. 

A  schematic  for  the  regenerative  heat  exchanger  is  shown  in  Figure  4. 

1st  Law: 


m  xhx(Px,  v,)  -m2h2(P2,  v2)  +  m3/i2(/?3,  vt)  -  rhXiP 4,  v4)  =  0 


thermal  characteristic: 


(7'I-7’4)-(r2-r3) 

mxhx  —  w  2/i2  —  UrAr  Ini  — — — —  — — — — —  =0 

.  i  22  r  r  i  (r,-r4)/(r2-r3) 


pressure  characteristic: 


Px-P2- 


'16 fl\  2 

—  m,(v,  +  v2)  =  0 

yKD  )n 
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Fig  4.  :  Regenerative  Heat  Exchanger  Schematic 
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Fig  5. .  Temperature  profiles  for  which  the  rate  equation  is  invalid 
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(l6fl]  .2 

Pz  PA  3  m3(v3  +  v4)  — 0 

)v. 

continuity: 

mx-m2  =  0  m3-m4  =  0 

The  temperature  characteristic  is  a  form  of  the  rate  equation.  Stream  1-2  is  assumed  to  be 
the  warm  stream  and  stream  3-4  is  the  cold  stream.  It  is  a  valid  representation  of  the  tem¬ 
perature  profile  through  a  heat  exchanger  for  profiles  that  exhibit  a  continuous,  exponen¬ 
tial  temperature  difference  along  the  flow  path.  Figure  5  shows  temperature  profiles  for 
which  the  rate  equation  is  not  valid:  discontinuity  in  the  profile  or  invariant  profile.  In 
instances  where  a  phase  change  occurs,  two  separate  rate  equations  should  be  used,  one 
for  each  side  of  the  discontinuity.  For  an  invariant  profile,  the  rate  equation  is  valid  by 
taking  a  limit  when  ATa/ATb  goes  to  0/0.  "UA"  is  the  component  heat  transfer  coefficient 
and  heat  transfer  area. 

Two  pressure  relations  are  given  here  for  the  friction  loss  phenomena.  For  model 
development,  the  pressure  drop  is  simply  represented  by  a  flow  duct  with  characteristic 
hydrr'-'lic  diameter,  flow  path  and  friction  loss  factor,  which  should  be  known  apriori  for 
a  given  component.  This  is  a  very  generic  representation  of  pressure  drop  due  to  friction. 
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The  segregation  of  the  characteristic  equations  for  this  device  into  temperature  and 
pressure  relations  does  not  imply  the  independence  of  heat  and  momentum  transfer.  On 
the  contrary,  they  are  intimately  coupled,  which  would  be  reflected  by  their  respective 
values  for  "f”  and  "UA”. 

An  ambient  heat  exchanger  is  shown  schematically  in  Figure  6.  The  constitutive 
relations  for  an  ambient  heat  exchanger  are  very  similar  to  those  for  the  regenerative  heat 
exchanger  with  the  second  flow  stream  being  replaced  by  a  quantity  of  heat  being 
rejected  to/  received  from  the  environment. 

1st  Law: 


mxhx(Px,vx)  -m2h2(P2, v2)  =  Q 


temperature  characteristic: 


mxhx 


-m2h2-  U4A4ln 


A^A' 


{T„-T2)-{T„-Txj 
(T„-T2)/(T„-T,)  . 


=  0 


pressure  characteristic: 


( 


Px-P2- 


16 fl' 

7tD3, 


"h(v,  +  v2)  =  0 


12 


continuity: 


m  \  -  m2  =  0 
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The  relations  for  the  ambient  heat  exchanger  are  of  the  same  form  as  those  for  the  regen¬ 
erative  heat  exchanger.  A  subtle  difference  exists  in  that  the  ambient  environment  is  rep¬ 
resented  by  a  temperature  reservoir  T„.  At  steady  state  or  for  very  large  reservoirs  ( i.e. 
Reservoir »  TOsys.tm)-  is  invariant.  T_,/ and  UaAa  are  also  defined  during  initialization. 


The  last  device  which  is  commonly  found  in  cryogenic  systems  is  the  throttling  or 
expander  valve,  shown  schematically  in  Figure  7. 

1st  Law: 


m  xhx(Px,  v,)  -m2h2(P2,  v2)  =  0 


characteristic  relation: 


hx-h2~  0 


continuity: 


m,-m2  =  0 

The  principle  characteristic  of  the  expansion  valve  is  that  the  entering  and  exiting  enthal¬ 
pies  are  equal  (although  the  enthalpy  through  the  valve  is  not  constant).  Note  that  the 
pressure  drop  across  the  valve  is  not  specified  as  a  characteristic  relation.  In  order  for  the 
existence  of  this  valve  to  make  physical  sense  in  a  refrigeration  system,  it  must  operate  in 
conjunction  with  a  pressure  elevating  device  (compressor).  Placing  a  required  pressure 
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ratio  on  the  expansion  valve  would  place  a  redundant  overconstraint  on  the  system,  since 
the  compressor  already  sets  the  system’s  principal  operating  pressure  levels  (see  Section 
2.2:Compatibility). 

4.4  Interconnectivity  Relations 

The  interconnectivity  relations  are  mathematically  trivial  but  form  the  critical  link 
to  the  algorithm.  By  equilibrating  the  independent  thermodynamic  variables  between  the 
outlet  of  one  device  and  the  inlet  of  a  connected  device,  thermodynamic  matching  is 
imposed  on  the  system.  Now,  all  elements  behave  in  a  highly-coupled  interactive  manner. 

Interconnectivity  places  three  relations  at  each  component  interface:  pressure,  spe¬ 
cific  volume  and  mass  flux.  Equilibration  of  pressure  and  specific  volume  ensure  a  ther¬ 
modynamic  match  of  the  working  fluid  and  equilibration  of  mass  flux  closes  the 
continuity  loop  around  a  closed  system.  Typically,  interconnectivity  relations  will  have 
form 

P2-P2  =  0  V2  V3  m2-m3  =  0 

4.5  Summary 

The  device  models  presented  here  are  self-contained,  problem-independent  equation 
sets.  The  next  task  is  to  shape  the  constitutive  relations  so  that  they  may  be  compatible  for 
assembly  of  the  [Kl]  matrix  and  the  (R)  vector. 
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5  Topology 


5.1  Introduction 

Prior  to  discussing  the  assembly  of  the  [Kl]  matrix  and  { R }  vector,  it  is  useful  to 
describe  how  the  topology  is  developed,  since  [K‘]  and  {R}  are  assembled  using  the  sys¬ 
tem’s  topolgy.  The  topology  establishes  in  mathematical  terms  which  elements  are  con¬ 
nected  and  how  they  are  connected  in  an  engineering  system.  Different  types  of 
interactions  yield  different  topological  structures  or  networks. 

The  most  efficient  numerical  method  of  representing  the  system  topology  is  with 
the  topology  or  incidence  matrix.  Such  a  matrix  is  a  tertiary  state  (1,  0  or  -1)  matrix  that 
relates  the  connectivity  of  nodal  variables.  Topology  matrices  are  commonly  found  in 
discrete  network  systems  such  as  electrical  network  circuits  exchanging  current  or  struc¬ 
tural  networks  exchanging  displacements.  In  the  Control  Volume  methodology,  topology 
is  used  to  describe  how  elements  exchange  mass,  energy  and  entropy.  Hence,  there  is 
more  than  one  interaction  at  each  node.  Figure  8  shows  a  typical  structural  network.  The 
springs  represent  elements  that  store  energy;  the  nodes  represent  equipotential  points 
between  storage  elements.  The  topology  matrix  for  this  system  is 

’  1  - 1  0  0' 

0  1-10 
Top[iJ ]  =0  0  1  - 1 

-10  10 
.-1  0  0  1  . 


37 


1 


The  columns  represent  nodes  (interconnective  elements)  and  the  rows  represent  storage 
elements.  In  the  CV  algorithm,  topology  is  defined  on  two  levels:  thermal  and  mass 
exchange.  Thermal  topology  describes  elements  in  thermodynamic  communication:  stor¬ 
age  elements  exchanging  heat,  work,  fluid  energy  and  entropy.  The  mass  exchange  topol¬ 
ogy  describes  the  arrangement  of  elements  exchanging  mass  and  momentum. 

Why  is  this  done?  It  is  principally  a  matter  of  mathematical  convenience  for  the  assembly 
of  the  appropriate  residual  functions  and  [Kl]  functions.  The  thermal  topology  is  used  to 
assemble  the  1st  Law  and  temperature  characteristics,  while  the  mass  topology  matrix  is 
used  to  assemole  the  pressure  characteristic  and  mass  continuity  relations.  The  distinction 
arises  due  to  the  presence  of  counterflow  heat  exchangers  that  have  one  1st  Law  relation¬ 
ship  and  one  temperature  characteristic  but  two  pressure  loss  functions  and  two  mass  flow 
relations.  In  such  an  arrangement,  the  two  separate  flow  passages  in  a  counterflow  heat 
exchanger  each  constitute  a  separate  pressure  loss  mechanism.  But  only  one  energy  inter¬ 
action  is  present. 

The  computer  will  generate  these  matrices  by  scanning  screen.  However,  the  ini¬ 
tialization  routine  will  require  two  separate  topological  generations,  one  on  the  thermal 
interaction  level  and  one  on  the  mass  transfer  interaction  level. 
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As  an  example  of  the  application  of  the  control  volume  toplogy,  consider  the 
Claude  cycle  shown  in  Figure  9.  Five  storage  elements  and  six  interconnective  elements 
(nodes)  are  shown.  This  system  will  have  five  1st  Law  relations  and  thermal  interaction 
characteristics,  but  six  momentum  interactions  since  there  are  six  distinct  flow  passages. 
The  thermal  topology  matrix  is 
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The  mass  and  momentum  exchange  topology  matrix  is 

given  here  for  the  same 

Claude  cycle,  but  the  schematic  difference  is  shown  in 

Figure  10. 

The  distinct  flow  cir 

cuits  of  the  regenerative  heat  exchanger  are  now  modelled  as  separate  entities  such  that 
the  corresponding  topology  can  be  used  to  generate  the  pressure  characteristic  relations 
and  the  mass  flux  relations. 


0  0  0  0  1 
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1-1900 
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The  convention  used  in  both  topological  representations  is  taken  as  flow  into  an  element 
as  positive  and  flow  out  of  an  element  as  negative. 
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6  Generation  of  the  Residual  Relations 


6.1  Introduction 

Chapter’s  4  (elemental  models)  and  5  (topology)  laid  the  mathematical  foundation 
for  the  constitutive  relations  that  serve  as  the  basis  for  generating  the  [K1]  matrix  and  the 
{R}  vector.  The  algebraic  equation  sets  presented  in  Chapter  4  must  be  manipulated  to 
conform  to  the  numerical  requirements  of  [Kl(xt)]{Ax}=(R(X;)).  The  [K1]  matiix  is  gener¬ 
ated  from  the  elemental  constitutive  relations  in  two  steps.  Since  [K1]  is  a  Jacobian  of  the 
system  residual  relations,  the  first  derivatives  (with  respect  to  the  global  variables)  of  all 
the  elemental  constitutive  relations  must  be  computed  and  these  first  derivatives  are  the 
constitutive  relations  of  the  elements  assigned  to  their  proper  [K‘l  locations.  This  is  done 
using  the  system  topology  information. 

6.2  Assembly  of  the  (R(x)}  Vector 

{ R }  is  the  array  used  to  determine  the  steady  state  operating  conditions  of  the  sys¬ 
tem.  It  is  also  assembled  using  the  system  topology  information. 

The  assembly  of  the  { R }  vector  is  much  simpler  than  the  assembly  of  the  [K'] 
matrix  (Section  6.4).  A  subtle  differences  between  the  generation  of  { R }  and  fK1]  is  that 
the  residual  relations  (R(x)}  are  multiple  function  equations(e.g.  a  1st  Law  has  several 
arguments  summed  together).  Hence,  their  recursive  relations  shall  contain  summa¬ 
tion  operators  to  cycle  through  the  nodal  variables  and  use  the  topology  to  filter  out 
entries  that  are  not  pertinent  to  the  specific  residual.  On  the  other  hand,  the  elements  in 

dAj 

the  [Kt]  matrix  are  usually  single-argument  entities  (e.g.^-)-no  summation  over  the  range 
of  nodal  variables  is  made. 
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As  an  example,  consider  the  1st  Law  relation  for  the  regenerative  heat  exchanger 


shown  in  Fig.  10: 


m1hl  -  m2h2  +  mAhA  -  m5h5  =  R 


The  thermal  topology  array  (row)  for  the  element  in  question  will  contain  (+1)  for  the 
flows  (from  nodes)  entering  the  element,  (-1)  for  flows  (to  nodes)  exiting  elements  and 
(0)  for  nodes  not  connected  to  the  element  under  consideration.  The  regenerator  in  Fig.  10 
corresponds  to  storage  element  3,  which  corresponds  to  row  3  in  the  thermal  topology 
matrix  Topr[3,j]  shown  on  page  42.  Thus, 

Topr[3,vl  =  [1  -10  1-1  0] 

R=  1  mJh1  Topr[3,y] 

i  =  i 

where  6  equals  the  number  of  nodes  in  the  system.  Expanding  the  summation  through  the 
appropriate  values  of  ’j’  (and  hence  TopT[3,j])  yields 

m,/t,  -rh2h2  +  mAhA -rh5h5  =  R 

Similar  combinatorial  operations  can  be  performed  to  generate  all  the  residual  functions. 
However,  not  all  constitutive  relations  have  such  well  regimented  structures  as  the  addi¬ 
tive  1st  Law.  In  such  cases,  the  numerical  value  of  the  topological  arguments  are  used  to 
modify  the  variable  indexes  and  parameters  where  appropriate.  This  is  apparent  in  the 
pressure  characteristic  constitutive  relation  (Section  4.3)  where  one  pressure  is  multiplied 
by  the  pressure  ratio  and  the  other  pressure  is  not. 
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The  topolgy  value  (1  or  -1)  can  be  used  to  ensure  that  the  pressure  ratio  multiplier  appears 
only  with  the  outlet  pressure.  The  details  of  the  recursion  relations  used  to  develop  the 
{ R }  vector  are  left  o  Appendix  1 . 

6.3  [K‘]  Matrix  Functions 

The  analytical  functions  that  form  the  entries  in  [K‘]  matrix  are  the  first  derivatives 
(with  respect  to  the  global  independent  variables)  of  the  residual  functions: 


K 


dl± 

dXj 


Taken  literally,  the  residual  functions  R,  must  be  known  prior  to  computing  K‘ir  This 
would  be  the  case  for  analytical  generation  of  [K1],  but  the  control  volume  method  is 
numerically  based  for  implementation  on  a  computer.  It  would  not  be  the  best  employ¬ 
ment  of  a  computer’s  potential  to  generate  and  store  all  the  system’s  residual  relations 
and  then  compute  their  derivatives.  It  would  be  prudent  to  store  the  unspecified  first 
derivative  constitutive  relations  (which  are  standardized  for  each  element  type)  and  then 
form  system  specific  constitutive  relations  from  these  standardized  relations  by  combina¬ 
tion  with  the  system  topologies.  This  is  precisely  the  what  is  done  here. 

As  discussed  in  Chapter  4,  all  the  possible  "R’s"  are  limited  to  1st  Law,  thermal  and 
pressure  characteristics  and  mass  continuity  relations  for  a  limited  collection  of  elements 
(Sect.  4.1).  X|  are  the  system  global  parameters,  namely  pressure,  specific  volume,  mass 
and  flow  rate  at  each  interconnection.  As  an  example,  a  1st  Law  relation  for  a  regenera¬ 
tive  heat  exchanger  is  developed  as  follows: 


where  p  represents  the  number  of  fluid  streams  entering  or  leaving  the  control  volume  for 
device  "i".  The  partials  of  R  with  respect  to  one  set  of  system  global  parameters  are: 


BRi  dR,  dR, 

dP’  dv/  drhj 


These  derivative  functions  have  the  following  standard  form: 


dRi  p  d/i: 

ap~=,?,m'ap‘]v 

dR,  p  dh. 


'dv, 


dR,  p 

~=  I  A,- 

om,  7=1  ’ 


where  the  right  hand  side  of  the  equations  contain  only  constitutive  relations  and  global 
variables.  The  use  of  the  topology  matrices  for  the  generation  of  the  [Kl]  terms  is  pres¬ 
ented  in  Appendix  1. 

As  mentioned  in  Section  2.3.2,  1st  Law  relations  with  heat  and  work  transfer  are  not 
incorporated  into  the  K‘  framework,  so  no  derivatives  need  be  computed  for  W  or  Q 
terms.  W  And  Q  computation  is  assigned  to  post-processing. 

Finally  a  numerical  (instead  of  an  analytical)  derivative  is  used  to  represent  the  line¬ 
arization  of  the  rate  equation  (heat  exchanger). 


6.4  Assembly  of  the  [K1]  Matrix 

6.4.1  Interaction  Regions  of  the  [K‘]  Matrix 

The  constitutive  relations  that  form  the  [Kl]  matrix  described  in  Section  6.3  are 
the  basis  for  the  formulation  of  the  [K‘]  terms. 

To  develop  [K‘],  it  is  necessary  to  visualize  the  matrix  consisting  of  separate 
interaction  regions  each  corresponding  to  a  type  of  interaction  as  described  in  Section 

2.3.2  and  Chapter  4:  1st  Law,  thermal  characteristics,  pressure  characteristics,  and 
mass  continuity. 


[K‘ (x)]  {x}  = 


1st  Law 

temp,  characteristic 
Pres,  characteristic 
.  mass  continuity  . 


{x} 


This  administrative  subdivision  of  the  [K‘]  matrix  is  necessary  due  to  the  two  differing 
types  of  interaction  that  lead  to  two  distinct  topology  matrices:  thermal  and  mass 
exchange.  The  vector  { x }  is  the  collection  of  the  independent  variables  relevant  to  the 
problem:  {P\,P2,  ...,P„,v„v2,  ...,v„,m„m2,  The  derivative  functions  pres¬ 

ented  in  Appendix  1  must  now  be  uniquely  assigned  to  their  appropriate  location  in 
[Kl]  using  the  topological  information. 

Since  there  will  only  be  three  types  of  independent  variables,  each  interaction 
region  within  the  [K']  matrix  may  be  further  subdivided  (in  an  administrative  sense) 
into  smaller  2nd-order  arrays. 
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Each  array  denoted  by  [  ]  represents  the  collection  of  1st  derivatives  of  a  constitutive 
relation  with  respect  to  one  of  the  independent  variables.  The  six  arrays  corresponding 
to  1st  Law  and  thermal  characteristic  are  identical  in  size  to  the  thermal  topology 
matrix,  used  to  assemble  the  arrays  from  the  constitutive  relations  of  the  elements. 
Similarly,  the  six  arrays  corresponding  to  pressure  characteristic  and  mass  continuity 
use  the  mass  topolgy  matrix  for  function  generation.  All  arrays  have  the  same  number 
of  columns  (since  each  column  represents  a  node),  which  ensures  the  columns  pro¬ 
perly  align  when  the  total  [K‘]  matrix  is  assembled. 

The  justification  for  this  arrangement  within  the  [K‘]  matrix  is  that  it  makes 
assembly  of  the  matrix  simpler.  Once  the  respective  toplogy  matrices  are  established. 
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each  storage  element  is  queried  for  its  applicable  constitutive  relations  (which  form 
{ R } )  and  [K‘j.  These  relations  are  combined  with  the  topology  through  recursive  rela¬ 
tions  to  define,  in  this  case,  the  appropriate  [K‘]  entry. 

6.4.2  Array  Assignments 

Each  array  shall  be  labelled  using  the  following  convention  for  ease  in  identify¬ 
ing  what  interactions  it  represents,  where  it  should  be  assigned  in  [K‘],  and  for  esta¬ 
blishing  the  recursions  used  to  produce  the  [K‘]  entires.  The  first  superscript  "t" 
indicates  (as  before)  that  the  functions  to  be  developed  consist  of  terms  of  the  system 
tangent  stiffness  matrix.  The  second  superscript  indicates  the  interaction  type  to  which 
this  array  refers(lst  Law  =1,  therm.  char.=T,  pres.  char=P,  mass  cont.=m.  The  first 
subscripts  are  the  array  location  indices.  The  outer  subscript  indicates  the  respective 
variable  of  derivation  (  P,  v  or  m).  As  an  example,  consider 

(*34  )\ 

Superscript  "1"  means  that  this  array  refers  to  a  1st  Law  relation.  Subscript  "P"  indi¬ 
cates  that  the  functions  inside  are  the  partial  derivatives  with  respect  to  the  indepen¬ 
dent  pressure  variable.  The  indices  [3,4]  indicate  that  this  particular  entry  corresponds 
to  the  third  1st  Law  relation  and  the  fourth  global  pressure.  This  labelling  has  been 
developed  so  that  each  inner  matrix  can  be  developed  independent  of  the  next  one 
using  the  recursive  relations  that  are  to  be  presented  in  the  following  paragraphs. 

Using  this  identification  scheme,  The  [K‘]  matrix  has  as  arrays 
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m\ 

ml 

ml 

ml 

ml 

ml 

ml 

_m: 

m: 

mi 

The  partial  derivative  functions  of  the  constitutive  relations  of  Appendix  1  (dis¬ 
cussed  in  Section  6,3)  are  combined  with  the  respective  topology  matrices  to  form  the 
inner  matrices,  which  form  the  global  [K‘]  matrix.  Each  row  within  an  inner  matrix 
corresponds  to  the  derivative  of  a  storage  element’s  constitutive  relations. 

As  an  example,  the  recursive  relations  to  generate  the  inner  matrices  for  the  1  st 
law  interactions  are: 


(Af')‘ =my^lvTopr[i,y] 

(Kti  =  rhJ^]PTopT[i>j\ 

(/^  =  /i/Topr[t,y] 


The  recursive  relations  for  all  inner  matrices  are  contained  in  Appendix  3. 

6.4.3  Summary 

The  [Kl]  matrix  has  been  subdivided  for  assembly  purposes  into  interaction 
regions  and  inner  arrays  corresponding  to  the  global  variable.  Each  inner  array  is  gen¬ 
erated  by  a  recursion  relation  (App.  2)  standardized  for  each  type  of  control  volume. 
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Any  inner  matrix  may  be  identically  zero  if  the  constitutive  relations  (for  that  array) 
are  not  functions  of  the  variable  used  to  compute  the  partial  derivatives.  Each  [K'] 
term  is  a  function  of  constitutive  relations  of  the  elements. 

Thermal  toplogy  is  used  to  generate  1st  Law  and  thermal  characteristic  functions 
while  the  Mass  topology  is  used  to  generate  the  pressure  characteristic  and  mass  conti¬ 
nuity  functions.  If  a  system  is  composed  entirely  of  two-port  elements,  then  the  mass 
and  thermal  topology  matrices  will  be  identical.  In  the  general  case  however,  the 
index  ranges  for  the  two  matrices  will  be  different. 
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7  Chronology  of  CV  Method 


7.1  Introduction 

The  formulation  and  solution  procedure  has  three  steps:  preprocessing,  generation 
(system  assembly)  and  reduction. 

7.2  Preprocessing 

The  first  step  in  the  preprocessing  operation  is  to  establish  the  control  volumes  of 
interest  in  a  closed  loop  system.  During  the  initialization  of  each  individual  storage  ele¬ 
ment  local  variables  are  assigned  based  on  the  type  of  element,  i.e.  a  standard  two-port 
element  (expander)  will  have  inlet  variables  Px  ,  v1  ,m'  and  outlet  variables 
P 2  ,  v2  ,m2.  A  four  port  element  will  have  Px,  ...,P4  etc.  The  local  variable  assignment 

sets  the  stage  for  the  local  to  global  variable  transformation  (i.e.  system  specification)  as 
discussed  in  Sections  2.3.1,  2.3.3  and  4.4. 

These  transformations  are  conducted  for  all  the  independent  variables  (P  and  v  as 
well)  at  the  interconnecting  nodes  only.  This  local  to  global  transform  has  two  functions: 
(1)  it  reduces  the  number  of  independent  state  variables  by  a  factor  of  two  thereby  making 
the  system  algebraically  unique  and  (2)  it  impresses  on  the  system  the  interactive  require¬ 
ments  of  the  individual  elements. 

When  the  elements  that  constitute  the  system  are  selected  and  connected  per  desig¬ 
ner’s  requirements,  the  interconnectivity  is  then  known  and  the  topology  can  be  estab¬ 
lished.  The  topology  is  then  used  to  generate  the  applicable  matrices  and  vectors. 
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7.3  System  Assembly 

The  assembly  revolves  around  the  individual  control  volume  as  a  basic  construc¬ 
tion  element.  Once  the  algorithm  determines  the  number  of  storage  elements  and  connec¬ 
tive  nodes  (via  the  topology  matrices),  it  will  establish  space  (memory)  sufficient  to  carry 
the  requisite  number  of  constitutive  relations  based  on  the  characteristics  of  each  control 
volume.  The  details  of  the  assembly  operation  were  discussed  in  Chapter  6. 

7.4  Reduction 

Reduction  is  the  actual  numerical  solution  of  the  [K‘(xk)]{Axk}  =  {R(xk)}=0  system. 
It  uses  a  simple  Newton-Raphson  method  to  successively  iterate  toward  the  equilibrium 
state,  {x^}  (the  index  "k”  represents  the  "kth"  numerical  value  based  on  thermodynamic 
state  { xk } ).  This  method  is  a  common  numerical  exercise  and  is  discussed  in  detail  in 
Appendix  4. 
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8  Conclusions 


This  effort  was  the  first  step  in  an  effort  that  is  to  develop  new  computer-based  tech¬ 
niques  for  the  simulation  of  cryogenic  engineering  plant  performance.  This  work  centered  on 
cycle  analysis  of  steady-flow  cryogenic  systems.  Simulation  is  the  extension  of  this  analysis 
to  more  generalized  operating  scenarios.  Analysis  is  indeed  difficult,  yet  it  is  only  the  first 
step  toward  simulation.  In  analysis  however,  some  important  foundations  can  be  established 
(or  alternatively,  bad  foundations  can  be  eliminated)  which  may  be  incorporated  into  simu¬ 
lation  models.  Such  matters  concern  variable  management,  idealization  of  the  physical  prob¬ 
lem,  mathematical  formulation  of  the  idealization  and  a  solution  algorithm.  There  are  some 
fundamental  obstacles  which  stand  in  the  way  of  achieving  a  simulation  model  for  a  thermal 
power  system.  These  obstacles  arose  during  this  effort  conspicuously  and  inconspicuously- 
certain  obstacles  remain,  some  were  overcome.  The  remarks  that  follow  address  in  retrospect 
the  extension  of  this  work  to  future  development  and  the  problems  in  simulation  in  general. 

Most  of  the  work  that  has  been  done  in  the  simulation  of  thermal  power  systems  has 
been  done  in  the  nuclear  power  plant  field  due  to  the  critical  need  for  understanding  of  plant 
behavior  before  embarking  on  construction.  However,  most  simulation  models  are  limited  to 
a  very  select  set  of  flow  circuits.  A  methodology  to  perform  analysis  (or  simulation)  of  arbi¬ 
trary  thermal  circuits  has  not  been  rigorously-developed  or  disseminated.  This  paper  pres¬ 
ented  a  method  of  analysis  for  arbitrary  thermal  power  systems. 

Another  fundamental  obstacle  in  analysis  and  simulation  of  power  systems  lies  in  the 
historical  practice  of  thermodynamics.  The  industrial  age  has  seen  the  expansion  of  power 
sources  (electrical,  automotive)  to  all  comers  of  the  earth.  Accordingly,  much  of  the  thermal 
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power  engineering  has  been  design  oriented.  As  a  result,  thermodynamic  models  for  pro¬ 
cesses  have  been  tooled  to  yield  design  data  based  on  some  required  performance  specifica¬ 
tion:  1st  Law  results  scale  sizes  of  systems;  efficiency  definitions  dictate  what  efficiency  is 
required  to  achieve  some  net  work;  rate  equations  define  area  and  heat  transfer  requirements 
for  a  specified  heat  transfer.  Such  models  form  the  cornerstone  of  time-proven  design  aids. 
Design  is  the  antithesis  of  analysis  and  simulation,  and  herein  lies  the  difficulty  for  the  analy¬ 
sis  problem:  a  plethora  of  good  design  models,  yet  a  dearth  of  useful  analysis  models.  What 
is  frequently  done  is  the  design  models  are  "reverse-engineered":  performance  is  determined 
based  on  specified  component  design  data.  The  practical  application  of  such  an  approach  has 
shortcomings. 

Design  models  are  basically  lumped-parameter  relations  (whereas  steady  flow  devices 
usually  exhibit  gradients.  These  lumped-parameter  models  are  useful  due  to  their  simplicity 
which  was  one  reason  for  implementing  them  in  this  method:  use  simple  models  to  build 
complex  systems.  More  complex  elements  will  be  handled  in  the  same  way  by  the  processing 
and  solution  algorithm.  However,  as  was  seen  in  Chapter  6  (and  in  the  appendices),  even 
these  simple  models  required  significant  algebraic  and  numerical  manipulation  to  fit  the  into 
framework.  Another  shortcoming  of  the  lumped-parameter  formulation  was  that  it  restricted 
the  class  of  problems  to  be  studied  (steady  flow  systems).  To  analyze  off-design  and  transient 
behavior  would  require  models  of  drastically  greater  complexity  which  are  not  expressible  as 
lumped-parameter  elements,  whose  behavior  may  or  may  not  be  expressible  in  closed  form 
analytical  functions.  In  effect,  the  analysis/simulation  problem  is  an  enigmatic  trade-off 
between  elemental  model  simplicity  and  algorithm  complexity:  a  numerically  based  analysis 
algorithm  begs  to  be  discredited  and  linear  while  thermal  systems  are  unquestionably  contin¬ 
uous  and  inevitably  non-linear.  If  element  models  and  their  interactive  behavior  could  be 
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linearized  and  discredited,  then  all  processes  could  be  simply  reduced  to  differences  amongst 
linear  functions  of  state  variables.  However,  it  is  well  known  from  the  study  of  thermody¬ 
namics  that  most  real  processes  are  path-dependent,  rendering  such  an  approach  ineffective. 

Future  study  may  investigate  similar  problems  in  different  fields.  For  example,  the 
finite  element  analysis  of  plastic  flow  and  failure  of  a  solid  is  a  problem  in  continuous,  non¬ 
linear  behavior  that  has  been  discredited.  However,  it  differs  from  the  thermal  power  prob¬ 
lem  in  that  a  power  system  may  have  a  various  type  of  constitutive  components  while  a  steel 
truss  is  constitutively  homogeneous  throughout  its  system. 

Another  approach,  as  a  corollary  to  the  CV  method,  may  be  to  impose  equilibrium  at 
the  elements  and  interpret  in  some  way  the  mismatch  at  connecting  nodes  (the  CV  approach 
imposed  equilibrium  at  the  nodes  and  allowed  the  elements  to  achieve  a  local  equilibrium).  It 
has  been  suggested  that  an  initial-value  approach  be  investigated.. .the  number  of  potential 
paths  is  still  unknown. 

This  author  wishes  to  illustrate  that  this  is  one  approach  to  the  problem,  an  approach 
which  sought  a  delicate  balance  between  element  simplification  and  algorithmic  complexity. 
At  the  expense  of  sounding  reflective,  the  author  feels  compelled  to  share  these  uncertainties, 
questions  and  apprehensions  for  the  edification  of  future  investigators  to  keep  this  line  of 
research  viable. 
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Appendix  1  Recursion  Relations  for  Generating  the  {R}  Vector 

As  in  the  assembly  of  [K‘],  two  topology  matrices  are  used  to  generate  {R ) :  thermal 
topology  to  generate  1st  Law  and  temperature  constitutive  terms,  and  mass  topology  to  gen¬ 
erate  the  pressure  characteristics  and  mass  continuity. 

1st  Law: 


R,  =  I  m  h  Topr[i,f] 

j  =  i 

The  upper  summation  limit  ”n"  corresponds  to  the  number  of  nodes  in  the  system. 
Thermal  characteristic  constitutive  relation: 

For  the  expander: 


fl,=  I  {hj  -T),/j  }Topr[i,/[ 

7  =  1 


'1' 


R2  =  1  {hj-T],hJs}TopT[i,j]  -  {Topr[t,/]  +  1} 


7  =  1 


\2J 


n 

^3  =  ^ { S(j + TopTii ,yi )s  -  S/} Topr[i,y] 

For  the  compressor,  the  usual  substitutions  apply. 
Pressure  characteristic  constitutive  relation: 

For  compressors: 


{Topr[j,y]  +  1} 
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(expander) 


*,■=  t[p  -iT^-i]Top„[i,y1 

j  a  J  c  -1 


For  fluid  flow  with  momentum  loss  (heat  exchangers): 


R,=  1 

j  =  i 


^  32/7^ 


Top, 


[/,/] 


Continuity  constitutive  relation: 


Appendix  2  Terms  of  the  [K‘]  matrix 

The  algebraic  system  described  in  Chapter  4  must  be  manipulated  to  conform  to  the 
numerical  requirements  of  the  [K‘]{Ax}={R}  solution  methodology. 

Partial  derivatives  of  1st  Law  Relations 

As  mentioned  in  Section  5.2,  1st  law  relations  that  have  environmental  interactions  are 
not  incorporated  in  K‘.  Thus,  a  first  law  relation  for  an  adiabatic  device  with  no  work  transfer 
can  be  written 

R.iXj)  =  I  m,hj  =  0 

j  =  i 

where  p  represents  the  number  of  fluid  streams  entering  or  leaving  the  control  volume  for 
device  "i".  P,  v  and  m  at  each  interconnection  are  the  independent  global  variables.  Hence, 
for  each  residual  relation,  the  following  partial  derivatives  are  of  interest: 


dR,  dR,  dR, 

dP/  dv/  drhj 

For  the  1st  Law  residual,  the  previous  equations  become: 


dR,  p  .  dhj 
dP~^midP] 


dR, 

dv, 


I 


i  =  i 


dR, 

drhj 
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These  recursive  relation  will  be  combined  with  the  topology  matrices  in  Chapter  7  to  yield 
the  actual  K‘  entries. 

Partial  Derivatives  of  Characteristic  Relations 

Expansion  engine  characteristic  functions: 

(hl-h2)-r\t(hu-h2s)  =  0 
h\~hu  =  0 

Direct  differentiation  with  respect  to  the  independent  variables  leads  to  the  linearized  charac¬ 
teristic  functions: 


dR{ 


dhy  dhu  dh^  ds^ 

"]v>  _T1,  dFx]v' +T1, 


&V 


ORy  dh2 

W2  =  ~W}V2 


dRy  _dhy  dhy  dh ^  ds * 

dVy  ~dVylP 1  ^'dVy  1  P  ^  ^  dsj  ^  dVy' ^ 


dR  y  dh2 

dv2  dv2]p* 

dR2  dhy  dhys 
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r 


II 

n 

bp2 

i 

cn  | 

£ 
<T>  \ 

3*2,. 

dp,' 

=  dF> 

M 
m  | 

1 

ds^ 

^2 


ivi 


3V]  dv,  '  dv[ 


dP2  dv2 


The  linearized  constitutive  relations  for  the  compressor  are  similar  to  those  for  the  expander 
engine  and  an  additional  characteristic  equation  for  the  pressure  elevation  term: 

Ra  =  P  ,  —  rEP2  =  0 

dR,_  a/?4_ 

W2~~Ke 


Regenerative  heat  exchanger  characteristic  equations: 


*i  =Pl~P2~ 


16/7 

TtD3 


"*i(vi  V2)  =  0 


12 
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'34 


R2  —  P 3  P4 


16/? 

tlO3 


m2Jv3  +  v4)  =  0 


/?3  =  rh{hx  -rh2h2  -  URAR  lnj 


(r,-74)-(72-7-3) 
{TX-T4)/(T2-T3)  . 


The  partial  derivatives  are: 


dRx 


dRx 

dP2 


dR ,  j 

(  32/? 

{ kD 3 

dRx 

f  i6 n 

3vi 

l^JtD3 

a/?, 

f  16/? 

dv2 

v7tD3 

dR2 

of-1 

br2 

r  32/7 

dm3 

^TtD3 

br2 

f  16/? 

3v4 

^TtD3 

\ 

m/Vi  +  Vj) 


\ 

rn]{vx  +  1) 

/ 


dR2 


=  1 


\ 

m3(v3  +  v4) 


\ 

m](l  +v3) 
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Special  attention  must  be  paid  to  the  linearization  of  the  rate  equation.  Obviously, 


dR3  dR2 

dm  i  1  dm2 


Partial  derivatives  with  respect  to  P  and  v  yield,  complex  functions  that  may  be  replaced  by  a 
numerical  derivative.  To  implement  this,  the  chain  rule  is  first  applied: 


/?3=^(^(p,v),r,(P,v),m,) 


^3  _  dg  dg  dT. 

dP~  hhidPi*  VTidPr' 


dR.^dgdh,  dg_d 7) 

dv,  dh,  dv,  p,  +  dTidvi  P' 


note  also  that 


dh, 


and 


dh, 

dv, 
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are  constitutive  relations  of  the  working  substance  and  would  exist  as  subroutines  external  to 
the  [K‘]  system,  evaluated  when  called  for  during  the  assembly  operation.  This  applies  as 
well  for 


and 


dIi 

dV; 


The  function  dg/dT j-  shall  be  evaluated  numerically  to  avoid  the  analytical  function  involving 

the  logarithm.  In  general,  the  numerical  derivative  for  a  function  of  multiple  variables  with 
respect  to  the  variable  X;  is 

dg(xux2, . . ..,xn)  g(xt,x2, . +  Ax,, . . .,.Q  - g (xhx2, . . . ..,xn) 
dx,  Ax, 


For  this  numerical  representation,  the  contribution  to  the  linearization  function  from  the  mihi 

terms  subtract  to  zero  since  h  is  not  an  explicit  function  of  T  in  this  postulation.  This  leaves 
only  the  non-linear  (logarithmic)  portion  of  the  rate  equation  to  be  evaluated,  denoted  f(T,): 


m=uAM 


<Ti-Tt)-{T2-T3j 
(T,-r4)/(r2-r3) . 


The  partial  derivative  of  the  rate  equation  residual  with  respect  to  (say)  P,  is 

8^_/(Tl)_/(r1(P,  +  AP„v1),r2,r3,r4)-/(r,(p1,v1),r;,r3,r4) 

ap,  ap,  a p, 


Other  partial  derivatives  are  straight  forward. 
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The  partial  derivative  relations  for  the  ambient  heat  exchangers  are  similar  to  the 
results  for  the  regenerative  heat  exchangers  with  the  expected  variable  and  parameter  substi¬ 


tutions. 


Expansion  valve  characteristic  equation: 


=  hl(Pl,vl)-h2(P2,v2) 


The  partial  derivatives  are  straight  forward: 


8/?, 

dPt 

3/?, 

dP2 


3/t, 

3/?, 

3v, 

8/2, 

~dv}]p' 

dh2 

dP2U' 

3/?, 

dh2 

dv2 

8v2  ^ 

Mass  continuity  relations 


Since  m,  are  independent  variables, 


dR  _  dR 
dP,~dv,=® 

dR  _  |  +  l(in-flow)  j 
dm,  l-l(out-flow)J 
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Appendix  3  Recursive  Relations  for  Generation  of  the  [K‘]  Matrix 

The  partial  derivatives  of  Section  6.2  are  combined  with  the  topology  matrices  to  form 
the  inner  matrices,  which  are  then  used  to  generate  the  global  [K‘]  matrix.  Each  row  within 
an  inner  matrix  corresponds  to  a  partial  derivative  of  a  constitutive  relation  for  a  storage  ele¬ 
ment. 


The  recursives  for  the  1st  law  interactions  are  : 


<=,3 


3P, 


Topf(i,/| 


(O'.  =  . 

m ,  - — 


dv, 


Topr[/J] 


(K‘)l  =  /tJop  T[i,j] 


The  recursive  relations  for  the  characteristic  equations  are  not  uniform  between  the  different 
devices.  For  the  work  transfer  devices. 


.T  _  r.  -dAy.  (dhjs.  ^;  +  Topr(,.y], 

(K]j)p  =  T°vT[i,f]{—]Vi-r\,{—]v  +■ 


ds 


BP, 


dsJ  + 


Topr[i,yl 


1  *  Top^li./I 


7+Topr(i,;li 

BP,  ' 


{Toprf/Jl  +  1}} 


Referring  to  Section  4.3,  the  recursion  for  the  compressor  engine  is  logically  extended  by 
replacing  -  1  j  with  (1  —  rjc) . 

For  the  expansion  valve,  the  temperature  characteristic  recursion  is 


66 


For  the  regenerative  heat  exchanger,  the  partial  derivative  of  the  temperature  characteristics 
become 


(K')^  =  A,Topr[tJ] 


(Ki  =  df(T,) 
dP, 


Topr[/,yj 


where  -zz-  is  the  numerical  derivative  defined  in  Appendix  1. 


-TT-  exists  as  a  subroutine 


with  standard  arguments  (T(Pj,  v,)).  When  the  topological  argument  is  non-zero,  the  corre¬ 


sponding  nodal  variable  will  be  passed  to  the  subroutine  computing  the  derivative  and  appro¬ 


priately  incremented  (AP„  or  Av,)  in  both  the  logarithmic  function  in  the  numerator  and  in  the 


denominator.  Similarly, 


(K'J  =  df(T,) 


3v, 


Top  T\iJ) 


Here  again,  the  thermal  topology  is  used. 

Pressure  characteristic  recursions  are  generated  by  combining  the  partial  derivatives  in 
Appendix  1  with  the  mass  topology.  The  pressure  relations  for  work  transfer  devices  require 
some  modification  to  the  constitutive  relation  to  make  it  compatible  as  a  numerically  gener¬ 
ated  partial  derivative.  For  pressure  elevating  devices,  the  pressure  elevation  is  always  speci¬ 
fied  on  the  outlet  term,  which  indicates  that  indicates  that  K'.j  must  be  able  to  distinguish  inlet 
from  outlet.  This  information  lies  in  the  topology  matrix  where  all  exiting  flows  have 
arithmatic  value  (-1).  To  implement  this,  the  pressure  residual  must  be  re-arranged  as  fol¬ 
lows: 


R^Pt-rS^O 


=  r?Pt-rlP. 


E‘  2 


4t°p„M  TTopj.,y] 

=  fc  P,-rF  I 


dRl 


|r;^f'-")Top  JU1  =(*:,), 
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and 


^l  =  Jr-i^.!«|Topjl,2l=(Ar;/ 


The  recursive  relation  for  a  pressure  characteristic  can  be  generalized  to 


<*;£= 


Top  m[i,j\ 


Application  to  compression  devices  is  simply  extended  by  replacing  rE  with  rc.  The  recur¬ 
sions  for  pressure  characteristics  experiencing  momentum  losses  are 

(K')Pp  =  Topn[i,j] 


f 


r/v;  +  v^i)  T°Pmt»j1 


(K‘)P  = 

'^v  II  J 


-16/7  i  2 

T  k(l+vy+TopMM)  TopJ/,/] 


The  mass  topology  matrix  is  used  to  modify  the  variable  index  on  the  specific  volume  so  that 
the  recursion  will  generate  the  exact  function.  Section  6.3  showed  that  when 


r  =  p,-p2- 


dR_ 


6-16/7^ 
V  nD'  ) 


yn\{\  +v,) 
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The  variable  index  "2",  which  is  the  linearizing  variable,  disappears,  but  the  variable  index 
”1”  remains.  The  matrix  variable  index  is  j=2,  so  j=2  must  be  modified  to  ensure  that  the 
variable  index  "1"  appears  in  the  function.  This  is  done  by  using  the  arithmatic  value  of  the 
mass  topology  matrix  to  mod  fy  the  variable  index  since  v2(an  outlet  variable)  has  arithmatic 
value  (-1). 

Finally,  the  mass  continuity  interaction  region  can  be  developed.  Clearly, 


The  mass  flux  partial  derivative  for  the  continuity  interaction  are  identically  the  mass  topol¬ 
ogy  matrix: 


(*')”  =  Topm[/J) 

where  "i"  corresponds  to  a  mass  continuity  constitutive  relation. 
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Appendix  4  Appendix  4:  Numerical  Solution  of  [K‘(Xj)]{Ax}={R(Xi)} 

Numerical  solution  (or  reduction)  of  the  [Kt(xik)]{Axk}=(R(xik)}=0  system  is  accom¬ 
plished  using  a  Newton-Raphson  scheme.  The  index  "k"  represents  the  "kth"  numerical  value 
based  on  thermodynamic  state  {xk}.  The  reduction  operation  consists  of  an  initial  estimate  to 
the  independent  variable  vector  { x° } .  { x° }  is  then  used  to  assign  a  numerical  value  to 
[K‘(x0)]).  The  product  [K^fAx0}  is  formed  which  is  the  first  estimate  for  the  system’s  con¬ 
stitutive  relations.  The  same  value  of  { x° )  is  used  to  compute  the  residual  relations  { R( x°) } , 
which  is  then  compared  to  {£},  the  error  vector.  Typically,  {R(x°)}>{e},  which  requires 
modifying  (iterating)  the  vector  {x}  until  { R(x°) } < { e } . 

The  (x)  vector  is  initialized  to  {x0}  by  the  user  based  on  some  educated  or  experience- 
based  guess.  When  { x°)  is  established,  [K‘(x)]  assumes  its  first  numerical  approximation, 
[K1,0].  Each  time  the  (R(xk)}  exceeds  the  error  vector  { e } ,  { x }  can  be  updated  by 

{xk+'}=-(sr(xtf' {#(**)} +  {**} 

{xk  +  '}={Axt}+{xij 

This  is  simply  Newton’s  Method  solved  for  (xk+I ),  the  improved  estimate. 

The  error  vector  { e }  is  used  to  check  for  convergence  toward  the  equilibrium  state 
{x^}  by  comparison  with  the  { R(xk) } .  Therefore,  { e }  should  reflect  the  scale  of  the  constitu¬ 
tive  relation  to  which  it  refers.  If  (e }  is  suitably  chosen,  it  will  ensure  quick  convergence 
with  acceptable  accuracy;  { e }  too  large  will  yield  rapid  convergence  but  poor  accuracy  while 


on  the  other  hand  {e}  too  small  will  lead  to  longer  convergence  times  while  generating  accu¬ 
racy  that  may  be  excessive  for  the  needs  of  the  analysis.  The  method  proven  useful  in  Dr. 
Russo’s  CAT  methodology  is  to  use  a  fractional  error  based  on  the  constitutive  relation  of 
interest: 


\R.(x°)\ 


where  "n"  represents  the  number  of  independent  variables.  For  a  small  system,  n  may  be  too 
small  to  render  desired  accuracy,  in  which  case  some  other  suitable  divisor  may  be  used.  This 
is  a  proven,  simple  and  consistent  system  scaling  technique. 

When  {x^}  is  attained,  the  post  processing  calculations  are  executed  since  all  the  inde¬ 
pendent  variables  are  known  and  all  the  post-processing  relations  are  functions  of  the  inde¬ 
pendent  variables.  As  a  minimum,  the  1st  Law  relations  for  non-zero  heat  and  work  transfer 
are  computed  to  complete  the  system  analysis.  Additional  computations  based  on  the  now 
determined  state  variables  may  now  be  conducted  as  specified  by  the  designer/operator. 
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